/*  scfcacos.cpp  by K.Tsuru  */
// function ID = 9115  since ver 2.18
/*****************************************************************
SComplex class
It returns arccos(z).
Cacos(z) = -i * Clog{z + Csqrt(z * z - 1)}.

Let z = x [+|-] i*y
arccos(z) = arccos(x [+|-] i*y) ([+|-]y >= 0)
= arccos[ (1/2){sqrt( (1+x)^2 + y^2 ) - sqrt( (1-x)^2 + y^2 )} ]
[-|+]i*arccosh[ (1/2){sqrt( (1+x)^2 + y^2 ) + sqrt( (1-x)^2 + y^2 )} ] (ref. Iwanami II p.224)
= arccos {(|1+z|-|1-z|)/2} [-|+] i*arccosh {(|1+z|+|1-z|)/2}

When y = 0 i.e. z = x, let arccos(x) = p, cos(p) = x then
|x| <= 1.0 must be satisfied.
******************************************************************/
#ifndef SN_H
#include "sn.h"
#endif

#define CacosXYMethod 0
#define CacosZ1Method 1

SComplex Cacos(const SComplex& z){
    if(Im(z) == 0.0){ // pure real z = x
        if(Dabs(Re(z)) > 1.0) {
           z.Real().SetError(z.Real().DOMAIN_ERR, "Cacos(z) with z = x(x > 1.0))", 9115);
        }
        return SComplex(Acos(Re(z)), 0.0);  // arccos(x)
    }
#if CacosXYMethod // 13.9sec
    SDouble x = z.Real(), y = z.Imag(), xp1sq, xm1sq, ysq(y*y), a, b, c, d;

    a = 1.0 + x;    xp1sq = a * a;  // (1+x)^2
    a = 1.0 - x;    xm1sq = a * a;  // (1-x)^2
    c = Sqrt(xp1sq + ysq);  // sqrt( (1+x)^2 + y^2 )
    d = Sqrt(xm1sq + ysq);  // sqrt( (1-x)^2 + y^2 )

    a = DsDiv(c - d, 2);
    a = Acos(a);

    b = DsDiv(c + d, 2);
    b = Acosh(b);   // b > 0
    if ( y.Sign(9115) > 0 ) b.ChangeSign(9115);

    return SComplex(a, b);

#elif CacosZ1Method // 13.87sec
    SComplex zp1(z + 1.0), zm1(z - 1.0);
    SDouble azp1, azm1, rp, ip;

    azp1 = Cabs(zp1);   // |1+z|
    azm1 = Cabs(zm1);   // |1-z|

    rp = DsDiv(azp1 - azm1, 2); // (|1+z|-|1-z|)/2
    rp = Acos(rp);

    ip = DsDiv(azp1 + azm1, 2); // (|1+z|+|1-z|)/2
    ip = Acosh(ip); // ip > 0
    if ( z.Imag().Sign(9115) > 0 ) ip.ChangeSign(9115);

    return SComplex(rp, ip);
#else // 16.3sec
    SComplex p = IU * Clog(IU * z + Csqrt(1.0 - z*z));

    p += MPi2();
    return p; // pi/2 + i*log{i*z+sqrt(1-z^2)}
#endif
}